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Abstract It is shown, by means of numerical simulations, that intercellular spiral waves of calcium can 
be initiated in a network of coupled cells as a result of a de-synchronization between Ca 2+ oscillations in 
two domains. No artificial heterogeneities need to be imposed to the system for spontaneous formation 
of spiral waves. The de-synchronization occurs near the interface of the stimulated region (which acts 
as a pacemaker) and propagates over the entire network. We also find the outcome of the collision of 
two spiral waves. 
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1 Introduction 

Oscillations and waves of cytosolic calcium (Ca 2+ ) have been reported in a large variety of cell types 
after stimulation by an extracellular agonist fTJ[2] . When stimulated, cells show a large variety of be- 
haviour which has been observed in their cytoplasm, including pulsating pattern, plane and spiral 
waves. Various experimental and theoretical works have been carried out to uncover the mechanism 
underlying these oscillations .ljSJIK]- in most models, Ca 2+ oscillations are the result of the auto- 
catalytic regulation by a phenomenon called calcium induced calcium release (CICR) by which Ca 2+ 
activates its own release from internal stores like endoplasmic reticulum. 

Spiral waves have been initially observed in the cytoplasm of cells such as cardiac myocytes [51 
[7], Xenopus Oocytes [8][9] and others [TOllTl"] . In a general manner, cytosolic spiral waves have been 
identified as a characteristic behaviour of excitable media [121113] ■ They can be initiated by spatial 
heterogeneity depending on the volume of the cell. That is in cardiac myocytes which are cells of 
small volume (1 mm in diameter vs. 100 fj,m length) the presence of an unexcitable region created by 
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nucleus [7] or by the existence of a region possessing a larger potentiality to release Ca 2+ in Xenopus 
Oocytes (very large cells with a diameter of up to 1000 /urn, much larger than myocytes |14j). 

An important characteristic of Ca 2+ oscillations is that Ca 2+ signals can propagate from cell 
to cell. They have been observed not only to propagate between cells of the same type (homotypic 
Ca 2+ waves) but also (at least in culture systems) between different cell types (heterotypic Ca 2+ 
waves). For instance, intercellular Ca 2+ waves have been observed propagating through ciliated tracheal 
epithelial cells [T51IT6] . rat brain glial cells [T7lfT8] and many other cell types [191120] . Furthermore 
intercellular Ca 2+ waves have been observed passing from glial cells to neurons and vice versa [211 
122] . For intercellular waves to propagate, two main types of coupling have been identified. The first 
mechanism involves the diffusion of a Ca 2+ mobilizing messenger through gap junction [23] and the 
second relies on paracrine signalling (24) involving the release of a messenger in the extracellular 
space, binding to receptors on the neighbouring cells, and activation of cytosolic Ca 2+ increase in 
target cells. A good review on calcium dynamics in general and on intracellular and intercellular 
waves in particular have been done recently by Dupont et al |25j . Particularly interesting are recent 
discoveries that intercellular waves are highly structured, with forms similar to the spatiotemporal 
organization of wave activity seen in Xenopus oocytes cytoplasm. In particular intercellular Ca 2+ 
waves in hippocampal slices have recently been observed to form spontaneous spiral waves [2l)ll271l2"8] . 
Numerical studies have been done to investigate how these spirals occur and how diverse parameter 
such as the intercellular conductance affects the process of spiral stability and breakup in an array of 
coupled excitable cells [28,29, 30 , 31 , 32 . It is now clear that intercellular Ca 2+ waves are a mechanism 
by which a group of cells can communicate with one another, and coordinate a multicellular response 
to a local event. Therefore, an understanding of the mechanisms underlying intercellular Ca 2+ waves 
is of importance, not only for a general understanding of intercellular communication, but also for 
the understanding of a wide range of specific processes such as mucociliary clearance, wound healing, 
mechanical transduction, cell growth, information processing and others. For example, in the wound 
healing response, if a monolayer of epithelial cells is mechanically damaged the resultant intercellular 
Ca 2+ wave sets up intercellular Ca 2+ gradients which influence the initiation and direction of cell 
migration [19] . 

In this paper, our goal is to propose another mechanism by which intercellular spiral waves can be 
initiated. We use an existing minimal model to numerically investigate the occurrence and propagation 
of intercellular spiral waves of calcium. The paper is organized in the following manner. Section [5] 
presents the two-dimensional array of coupled cells while in section [3] we give the results of the 
numerical simulations. Conclusion is given in Section 2] 



2 Model 

To model intercellular calcium waves, two aspects must be considered: intracellular dynamics of cal- 
cium and coupling between cells. Two types of theoretical models have been developed so far: the 
spatiotemporal models and the temporal models. The spatiotemporal models take into account the 
fact that intracellular calcium waves are spatially distributed in cells. However in some cell types, with 
small diameter, such as hepatocytes (10-20 \im) [331134] and pancreatic acinar cells(10-20 fim) [35 , 
in which the intracellular propagation velocity is at the order of lO/im.s -1 while intercellular speed 
is around 120 fim.s -1 [3S], the spatial intracellular aspect can be neglected. Thus, the dynamics of 
the cells can be approximated by a set of ordinary differential equations. An interesting study by 
Tsaneva-Atanasova et al. [35] investigated in pancreatic acinar cells the temporal and the spatiotem- 
poral models. Although the point-oscillator model can not explain all the phenomena exhibited in the 
cell networks such as synchrony, it has been found in Ref. [35 that it gives a reasonably accurate gen- 
eral picture when studying wave propagation. That is why we use the temporal model in this paper. 
However, to take into account the fact the signal wave takes a time to travel along a cell in some 
biological organs, we also carry out the numerical simulation of the model with a time delay added. 
This delay represents the time taken by a signal to propagate inside a cell before moving to the next 
cell. This constitutes a substitute or equivalent to the spatiotemporal model. 

The model and numerical computations are based on the minimal model of Dupont et al [J] based 
on CICR that was originally designed to model intracellular Ca 2+ oscillations in a cell and recently 
used to model intercellular propagation of Ca 2+ waves in a ID network of diffusively coupled cells [55] . 
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A two-dimensional network of coupled cells is considered in this paper. As in Ref. [36], we assume 
that cells are coupled together by a bidirectional paracrine coupling. Let us consider Xij as the Ca 2+ 
concentration in the cytosol and yij the Ca 2+ concentration in the internal store. Therefore, the cell 
defined by the i th and j th coordinates is described by the following set of equations (See Ref for 
detailed on the modeling): 



dx j j 
dt 



aij~V2,ij+V 3 ^, J +kfy l j-kx l j+f3 1 V 1 (x i+1 j-2x l:J +x i - 1 j)+f3 2 V 1 (xij +1 ~2x i:] +Xij- 1 ) (1) 

' / ' /; ' Vw-Vw-kfVij (2) 



dt 

with i = to N and j — to M. Node indicates a different cell 
In these equations, 

Vb + bV\ if excited 

Vq if not excited 



represents the term characterizing the excitation state of a cell. Vq represents a constant influx of Ca 2+ 
from the extracellular media to the cytosol whereas bV\ represents an external excitation which can 
be due to a hormonal stimulus binding to receptors in the extracellular membrane of the cell. The 
binding to the receptors causes opening of ionic channels. fi\ and fa represent the coupling constant 
respectively in the x and in the y direction. 

represents the speed of Ca 2+ pump from the cytosol to the internal store. 

v _ Vm3^fjj/fj 

represents the speed of Ca 2+ liberation from the internal stores to the cytosol. The activation of this 
process is provoked by the Ca 2+ itself characterizing the CICR process. The Ca 2+ extrusion from the 
cytosol to the extracellular media is taken into account by the term kXij. The Ca 2+ can also pass 
from the internal stores to the cytosol via the passive flux given by the expression kfyij. 

For the numerical simulation, no flux boundary conditions are used at the edge of the domain, 
defined as: 

Xo.j = XN+l,j = 
X t fl = X it M+l =0 

We solve the model equations in an array of 50 by 50 cells shown in Fig. [T] (extension to larger arrays is 
used when necessary); parameters chosen for the simulation are listed in Table [TJ they have been taken 
in Ref. [1] . The values of the coupling strength used in this work have been deduced from the study of 
Ref |3B] and in a recent study by Gracheva and Gunton [37]. Equations (1) and (2) are extended later 
(in the next section) to include time delay representing the duration of signal propagation in a cell. 



3 Results 

We perform the numerical simulation by integrating the model equations together with the boundary 
conditions on a grid of 50 by 50 cells shown in Fig.[T]. The fourth order Runge-Kutta algorithm is used 
with the time step equal to 0.001. For most results, a circular region with radius 5 cells is assumed to be 
stimulated at the center of the domain. This stimulated domain is defined as (i — io) 2 + (j — jO) 2 — 5 2 , 
where («o Jo) is the coordinate of the cell at the center of the excited domain. 

When the excitation act as a Dirac (for a localized mechanical or electrical excitation), the excited 
cells show a pick of calcium which can propagate to 2-4 neighbouring cells. However, when using a 
continued excitation to investigate spiral waves occurrence, one sees that when the degree of excitation 
is weak, a 1:1 locking is observed between cells of region 1 and region 2, therefore circular concentric 



4 



Wilfred D. Kepseu et al. 



Table 1 typical simulation constants for the minimal model 



Parameter 


Value 


k 


2s" 1 


kf 


1.0s" 1 


k> 


1.0/xM 


k a 


0.9/xM 




2.0fiM 


Vb 


1.3HM.S- 1 


v 1 


7.3/iM.s" 1 


V m 2 


Q~.0nM.s~ 1 


V m 3 


500.0fiM.s~ 1 


Pi 


0.50 


ft 


0.50 




Fig. 1 Typical geometry of system used to study intercellular spiral Ca 2+ wave. Only two spatial dimensions 
are considered. Outer square represents a 50 by 50 cells. First inner circle (region 1) is region in which cells are 
stimulated. Second inner circle (region 2) is a small region in the non excited zone choose in other to observe 
the behaviour of non excited zone. In that frame, region 1 is a circular region of radius 5 cells. 

waves are observed to propagate in the array. However, when the stimulation is strong, 1:1 locking is 
no more satisfied as seen in Fig. [3J where cells in the outer region of stimulation (region 2) exhibit 
weakly chaotic time behaviour, whereas cells in the stimulated region (region 1) show strong chaotic 
oscillations. Broad intercellular waterfronts with a width of 4-8 cells occur. These intercellular waves 
originate in the intermediate domain between the two regions and propagate in curvilinear and spiral 
patterns as shown in Fig. [3] Spiral waves are the result of waves de-synchronization between cells in 
the region 1 and region 2 which provoke a lateral instability at the intermediate area in between (which 
acts as a pacemaker). 

In a general manner, Ca 2+ spiral waves occurrence does not depend on the shape and dimension 
of the stimulated region (region 1). However, when the dimension of the stimulated area decrease, 
one need to increase the coupling between cells in order to observe spiral occurrence. The numerical 
simulations also show that varying the area of the stimulated region can induce a change in the rotating 
direction of the spiral waves as shown in Figure |H It is also observed that when the degree of excitation 
is further increased (b > 0.5 for the same choice of parameters), it is possible to observe at different 
times, calcium puffs propagation (Fig.[SJi), propagation of concentric waves (Fig. [5}d) and spiral waves 
(Fig. [Sp) which can suddenly occurs at a time further in the network. Calcium puffs are characterized 
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Fig. 2 Calcium oscillations of a cell in: (a) the stimulated region (region 1), (b) the outer region of stimulation 
(region 2) with b = 0.5. 



by small intercellular waves limited to 5-15 cells. The observation of the temporal behaviour of Ca 2+ 
concentrations of cells shows that there is still chaotic evolution in the outer region of excitation 
whereas cells in the stimulated region shows a burst (Fig. [5]) characterized by a higher pick of calcium 
and return to the equilibrium state. 

Extending the array dimensions does not have an effect on the appearance of spiral waves. Con- 
sidering for instance an array of 80x80 cells, when stimulating another region of radius 5 centered at 
the cell (i = 75, j = 25), one can see that when two curvilinear wavefronts collide, they annihilate at 
the point of contact. After collision, portions of the waves that have not collided merge and continue 
to propagate tangentially from the site of collision (Fig. [?])■ This behaviour is a standard property 
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X 

Fig. 3 Spiral Ca 2+ waves propagating in the network for b = 0.50. 




6 12 18 24 30 36 42 48 

X 



Fig. 4 Spiral waves showing the change in the rotating direction of the spiral. Figure obtained for the same 
choice of parameter when stimulating a region of radius r = 3 and 6 = 0.5. 

seen in excitable systems and which has already been reported experimentally in hippocampal slices 
cultures [53] and is a standard property of excitable media. 

Since the calcium wave in some cell types propagate slowly, one needs to treat the cell as a spatially 
extended system in which the signal propagates from one side to the other. In order to take into 
account this fact without using partial differential equations, we introduce a time delay representing 
the duration of propagation of calcium wave in each cell. By so doing, the set of equations (1) and (2) 
becomes (we consider here the case where the wave propagates in only one direction): 



= a>ij - V 2 ,i,j + V 3:i j + kfyij - kxij + PiVi(xi-ij(t - t) - Xij) + f3 2 V 1 (x ij - 1 (t - r) - xy) (4) 
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Fig. 5 Different patterns observed when the degree of excitation is high enough (b = 0.54). (a) Calcium puffs 
propagation in the network obtained at t = 17.19s. (b) Concentric waves propagation obtained at t — 39.09s. 
(c) Spiral waves obtained at t = 58.07s. 
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time(s) ( b > 

Fig. 6 Calcium oscillations of a cell, (a) In the stimulated region where cells show a burst, (b) In the outer 
region of stimulation, cells show a chaotic behaviour for b — 0.56. 



-J^ = V 2 ,i,j - \ - k f y i}j (5) 

with i = to N and j = to M. 

In this model, r is the time lapse necessary for the signal to propagate from one edge to the 
other of a cell. We consider in this case a unidirectional coupling. It is observed as shown in Fig. [5] 
that spiral still arise in the network, however, the de-synchronization occurred only in one side of the 
excited region. The effects of time delay on the speed of calcium wave have also been analyzed. To 
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(b) 

Fig. 7 (a) Collision of two spiral waves in the array, (b) Annihilation and propagation of the portion of the 
waves that have not collided for b = 0.5. 

determine the speed of calcium signal propagating in the array, we look at the wave propagation along 
one direction and we take the times t% and tk when the signal arrives at two different sites i and k and 
the speed is the quantity k — i divided by tk — U. Figure |H] shows that the speed of Ca 2+ oscillations 
in the network decreases when the time lapse r increases and increases with the coupling coefficient 
j3±. This is understandable, since it is known that increasing the time lapse r implies that the signal 
takes more time to propagate across a cell. 



4 Conclusion 

It is well known that a circular front that breaks in an asymmetric medium can initiate a spiral. 
This assumption has been used to show that intracellular spiral calcium waves can be initiated by 
an unexcitable region in cardiac myocytes^ or an excitable region in Xenopus oocytes[Tl]. Also, it is 
known that intercellular spiral calcium wave can be initiated by simulating the release and propagation 
of inositol 1,4,5-trisphosphate in several homes in Xenopus leavis[27] or by the presence of a bolus (in 
the shape of a line) of Ca 2+ which is placed directly behind the refractory region of a wave[5B]. In 
addition to other studies aimed to investigate the origin of spiral calcium waves, in the present work, 
it is shown that intercellular spiral calcium wave can also be simply initiated in a network of coupled 
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Fig. 8 Spiral waves observed when considering the delay between cells with b = 0.5. 




iH 

Fig. 9 Speed of propagation of Ca 2+ waves obtained when varying the time lapse /3i = 0.5, f$% — and 
6 = 0.5. 



cells as a result of the de-synchronization at the interface of an excited region and a non-excited region. 
Also, the outcome of the collision of two spiral waves has been found. This work complements a recent 
study using the spatiotemporal description of calcium flow and where it was found a mechanism of 
spiral generation at the interface between a pacemaker region and an outer region owing to the chaotic 
pulse transmission at the interface [55]. 
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